\section{Introduction} 

Deterministically modeled chemical reaction networks are of interest in the
field of systems biology and elsewhere. Models based on the assumption of mass-action
kinetics depend on the kinetic parameters of the system. However, measuring
these parameters empirically is difficult and error-prone.  Hence we focus on
establishing properties that are independent of the particular choice of
parameters.

In this work we address the issue of existence of positive equilibria i.e.\ 
concentrations of species that will stay constant under the system's dynamics.
The particular case we tackle is where the directed graph of chemical reactions
forms a strongly connected component and the reactions conserve mass.


\subsection{Background}
 A mathematical theory for this problem, also known as \textit{Chemical
 Reaction Network Theory} (CRNT), has its roots in the seminal work by Fritz
 Horn, Roy Jackson, and Martin Feinberg in \cite{GMAK, uniqueEPandLyapunov, necc-suff-CB, deficiency0, deficiency1}. We will
 build on the notation used in these works and summarized by Gunawardena et al.
 in \cite{gunawardena}. The system consists of a collection of species reacting
 collectively in some combination to give another combination of species in a
 network of chemical reactions. We will restrict our analysis to systems that conserve mass. Let $\mathcal{S}$ be a set of $m$ species and
 $\mathcal{C}$ be a set of $n$ complexes. The relation between species and
 complexes can be written as a non-negative matrix $Y\in\R^{m\times n}$, where column $y_j$ represents complex $j$ and 
 $Y_{i j}$ is the multiplicity of species $i$ in complex $j$. For example, the
 multiplicity of the species NaCl in the complex (H$_2$O + 2NaCl) is 2.  

A reaction network is represented by the underlying weighted directed graph
$G(V,E)$ with each node in $V$ denoting a complex, each directed edge
$i\rightarrow j$ denoting a reaction using $i$ to generate $j$, and the
positive weight $k_{i\rightarrow j}$ being the reaction rate.  The matrix $A
\in \R^{n \times n}$ is the weighted adjacency matrix of the graph, where
$A_{ij}=k_{i\rightarrow j}$.  Define $D := \diag(A\1)$, where $\1$ is the
vector of all ones, and $A_k := A^T-D$. The vector for concentration of species
in a reaction network is given by $c\in\R^m$. A nonlinear function
$\psi(c):\R^m_+\rightarrow \R^n$ capturing the mass-action kinetics is given by
\[\psi_j(c) = \prod_i\,c_i^{Y_{ij}},\] where the system of ordinary differential equations governing the
reaction network can be written as \[\dot{c} = YA_k\psi(c).\]
\noindent Hence, an equilibrium concentration for a chemical
reaction network is any non-negative vector $c^*\in\R^m$ such that
$YA_k\psi(c^*)=0$. Equivalently, a pair $(c^*,v^*)$ with $c^*\in\R^m$
and $v^*\in\R^n$ will be an equilibrium if it satisfies the conditions
\begin{align}
  YA_kv^* &=0 \label{fb}\tag{FB} \\ \psi(c^*) &= v^* 
  \label{mak}\tag{MA}
\end{align} or \emph{flux-balance} and \emph{mass-action},
respectively. We can find an equilibrium $c^*$ by finding a vector $v^*$
that satisfies $\eqref{fb}$ and also satisfies \eqref{mak} for some
vector $c^*\in\R^m$.

Observe that if $v^*$ is a positive equilibrium, then 
\begin{align}
  Y^T\log(c^*)= \log(\psi(c^*))= \log(v^*).
  \label{mak-alt}
\end{align} 
We refer to this alternative condition as the logarithmic form of the
mass-action condition \eqref{mak}.

In this notation, systems with the property of mass conservation can be characterized by the following definition. 
\begin{defn}
A chemical reaction network is $\emph{mass conserving}$ if and only if there exists a strictly positive vector $e\in\R^m$ such that 
\begin{align}
 e^TYA_k=0,
  \label{consis}
\end{align}
where $e$ denotes the molecular weights of the species (or atomic
weights if the species are elements).  
\end{defn}

We further restrict our analysis to a class of weakly reversible networks.
\begin{defn}
A chemical reaction network is \emph{weakly reversible} if for any complex $i \in\mathcal{C}$ there exists at least one reaction producing $i$ and one reaction consuming $i$.
\end{defn}
As suggested by simple examples like a single non-reversible reaction,
reversibility at least in a weak sense is a prerequisite for steady states with
strictly positive concentrations for all species.

The connectedness of the networks is captured in the following definition of a
terminal-linkage class.  
\begin{defn} A \emph{terminal-linkage class} is defined as a set of
  complexes $\mathcal{L}$ such that for any pair of complexes $(i,j) \in
  \mathcal{L}$ there exists a directed path in the graph $G$ that leads from
  $i$ to $j$.  
\end{defn} 
A reaction network that consists of exactly one terminal-linkage class is called a \emph{single terminal-linkage network}.

Another useful definition for a network is a stoichiometric subspace.
\begin{defn} A \emph{stoichiometric subspace} $S$ is the subspace defined by
  the span of vectors $y_{j}-y_{i}$, where $y_{j}, y_{i}$ are complex vectors for each reaction pair $i \rightarrow j$ in
  the network.  
\end{defn} 
In one of the early works, Horn and Jackson \cite{GMAK} analyzed
  mass-action kinetics and defined a class of equilibrium points, called
  \emph{complex-balanced equilibrium}, and systems admitting such equilibrium,
  \emph{complex-balanced systems}. These systems are shown to
  satisfy the \emph{quasi-thermostatic} and \emph{quasi-thermodynamic}
  conditions regardless of the kinetic rate constants. Following this, Horn
  \cite{necc-suff-CB} also proved necessary and sufficient conditions for
  existence of a \emph{complex-balanced} equilibrium. In
  \cite{uniqueEPandLyapunov}, Feinberg and Horn used the existence of a
  Lyapunov function to show the uniqueness of the strictly positive steady
  state in each stoichiometric compatibility class, which is equivalent to
  specifying all the conserved quantities of a system. Later Feinberg
  \cite{deficiency0, deficiency1} proved two theorems, now famously known as
  Deficiency 0-1 theorems, that provide the analysis of strictly positive
  steady states for a class of networks with deficiency $0$ or networks with
  deficiency $1$ but with each terminal-linkage class having deficiency less than $1$.
  The \emph{deficiency} of a network is defined as $\delta = n - t - s$, where
  $t$ is the number of terminal-linkage classes and $s$ is the dimension of the
  stoichiometric subspace, also known as stoichiometric compatibility class.
  For this restricted class of networks, the existence of a positive steady
  state is given by Perron-Frobenius theory for a strictly positive
  eigenvector. Other work in this area is from the perspective of dynamical systems and aimed toward proving two open conjectures:
  \emph{Global Attractor Conjecture} and \emph{Persistence Conjecture} \cite{Anderson-GAC}. Another approach using parametrized convex optimization to compute a non-equilibrium steady state is given in \cite{fleming-opt}.

To the best of our knowledge, the vast majority of CRNT research studies \textit{complex-balanced systems}, which, by definition, are assumed to admit a
\textit{complex-balanced equilibrium}. In the notation above,
complex-balanced equilibrium holds when the vector of concentrations satisfies
$A_k\psi(c)=0$, i.e., the vector $\psi(c)$ belongs to the null space of $A_k$.
It can be shown that a network with linearly independent complexes will have
deficiency $\delta = 0$. However, if some of the complexes are linearly
dependent (as shown in the example in Section 3), there are systems that are
not complex-balanced but admit equilibrium points with vectors $\psi(c)$ such
that $A_k\psi(c) \neq 0$ but $A_k\psi(c)$ belongs to the null space of $Y$. The
condition of \emph{complex-balance} is sufficient for \emph{thermodynamic}
consistency. However, as shown in \cite{GMAK}, it is not necessary. Also, for inhomogeneous systems with external force, the \emph{thermodynamic}
consistency conditions change and \emph{complex-balance} does not necessarily
hold.

\vspace{0.2in}
For the class of networks that admit thermodynamically
consistent equilibria for the homogeneous case, we would like to extend the results to the
inhomogeneous case, and for such networks, prove results on existence and
computability of equilibria that are not necessarily complex-balanced. In this
paper, we define a function based on a convex optimization problem.
The objective function of our formulation is very close to the Helmholtz
function defined in \cite{GMAK}. The fixed point of this function
gives the required equilibrium solution.  
We prove the existence of a strictly positive steady state for any weakly
reversible chemical reaction network with a \emph{single terminal-linkage class}. We
strongly believe that this can be extended to systems with \emph{multiple terminal-linkage classes}, as supported by our computational results for randomly generated networks. Section
3.2 gives a detailed analysis of a toy network to emphasize our extension.
